Self-propelled non-linearly diffusing particles. Aggregation and continuum description. 
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We introduce a model of self-propelled particles carrying out a Brownian motion with a diffusion 
coefficient which depends on the local density of particles within a certain finite radius. Numerical 
simulations show that in a range of parameters the long-time spatial distribution of particles is 
non- homogeneous, and clusters can be observed. A number density equation, which explains the 
emergence of the aggregates and indicates the values of the parameters for which they appear, is 
derived. Numerical results of this continuum density equation are also shown. 
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The tendency of living individuals to aggregate is an 
ubiquitous phenomenon in Nature that constitutes a cen- 
tral problem in different areas of the natural sciences. 
The usually discussed examples in the literature arise 
from very diverse contexts and comprise a wide range 
of spatiotemporal scales: human crowds, animal groups, 
cell populations or bacteria colonies. In this context dis- 
tinct situations are studied, including the dynamics of 
groups of individuals moving coherently resembling the 
behavior of fish schools, bird flocks or insect swarms , 
the pattern formation in populations of bacteria Q, or 
the patchy structure of plankton mediated by hydrody- 
namic driving Two levels of descriptions are usually 
performed, one considering the discrete particle dynam- 
ics, and the other accounting for the large spatiotemporal 
scales in terms of an evolution equation for the number 
density of particles. In the best of the cases, the last is 
derived from the more fundamental particle dynamics. 

From the theoretical point of view one of the main 
questions one can address is about what causes aggre- 
gation to form jij. In order to answer this, and start- 
ing from the discrete particle dynamics, many differ- 
ent mathematical models have been introduced which, 
mainly in the physics literature, are characterised by 
their simplicity but also because they retain the basic fea- 
tures underlying the clustering phenomenon. Essentially 
these models assume that i) individuals are self-propelled, 
and ii) they interact with the environment and/or with 
other particles through attractive and repulsive forces , 
or rather, by modifying its velocity because of commu- 
nications with neighbours [f|. Note, however, that even 
a simpler mechanism that cannot be classified within ii) 
and that gives rise to clustering has been introduced in 



ties, reducing, its random motility (diffusion coefficient) 
depending on the total number of particles in a neigh- 
bourhood of spatial range R. This response of the parti- 
cles to their local environment can be interpreted as an 
effective deceleration because of the presence (via colli- 
sions, interchange of chemicals, excluded volume effects, 
demographic pressure, etc..) of other particles. The in- 
teraction radius R appears frequently in the modelling of 
biological systems to take into account visual or hearing 
stimuli, chemical signaling, and other kinds of interac- 
tions at a distance. R is considered to be a fixed number 
in this work but in some biological species it depends 
on the environmental conditions. Just as an example, 
we mention the detection distance of preys for a type of 
zooplankton In our model, despite the particles are 
moving randomly and no attractive forces among them 
are considered, they tend to cluster. It will be obvious, 
when the coarse-grained number density equation of the 
model is presented, that this is because of the nonlin- 
ear character of the diffusion of the particles. Therefore, 
the aim of this work is to introduce this new type of 
mechanism, the local modification of the motility by the 
crowding of the surroundings, at the level of discrete par- 
ticle dynamics, and show under which conditions it gives 
rise to clustering. Then obtain its continuum description 
and interpret the model and its clustering instability in 
terms of this. The manuscript goes precisely along these 
lines: first we introduce the particle model and study 
it numerically, then we derive its continuum description 
and study it analytically and numerically. 



In this letter, following this line of searching simplic- 
ity, we introduce a basic model going through steps i) and 
ii) characterised because the aggregation phenomenon is 
due to a new type of mechanism, where no attractive nor 
repulsive forces among particles are taken into account. 
The model considers a fix population of particles mov- 
ing Brownianly. At every time step, any particle modi- 



Let us consider N pointwise particles initially dis- 
tributed randomly in a two-dimensional system of size 
L x L with periodic boundary conditions. At every 
discrete time step t the positions of all the particles, 
Xi(i) = (xi(t),yi(t)) (i = 1,...,JV), are updated syn- 
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chronously as follows: 

Xi(t + At) = Xi {t) + J ™^Lrg{t), 

y i{ t + At)= yi (t) + ^-^^r,V(t), (1) 

where Dq is a constant named in the following constant 
diffusivity or motility (this last name is the usual one 
for cells, bacteria or plankton organisms), At is the time 
step which we take equals to 1 in the numerical simula- 
tions, rj = {rj x ,7f) is a white noise with zero mean and 
correlations < T]f(t),r]j(t') >= SabSijStv- N R {i) denotes 
the total number of particles at distance less than R of 
particle i including itself, and p is a positive real num- 
ber. The model clearly states that particles are mov- 
ing Brownianly with diffusion coefficient inversely pro- 
portional to the number of neighbours within a range 
R, Di — D /(N R (i)) p , being Di the effective diffusiv- 
ity of particle i. In this way, if the neighbourhood of a 
particle is poorly crowded it diffuses quickly. The func- 
tional dependence of Di resembles the standard models 
of density-dependent dispersal of insects or animals 0, 
where the behavioral parameter p determines this depen- 
dence. Note, however, that in in these models, at differ- 
ence with the one introduced in this work, the diffusivity 
increases with the density of particles. A distinct type 
of decreasing function of Di with Nr(i) could be consid- 
ered, and will be discussed later on at the level of the 
density equation. 

First, we pay attention to the limits R — > and i? — > L 
(because of periodic boundary conditions the total sys- 
tem size limit for R is in fact smaller than L but we keep 
this notation just for clarity). In the first one, Nr(i) = 1 
for all i so that the particles diffuse with the constant 
motility, Di = Dq. In the second limit again all the par- 
ticles diffuse with the same diffusion coefficient given by 
Di = Dq/N, that is negligible for large N. Of course we 
do not expect any grouping behavior in these two limits 
and R will be considered in the rest of the paper to be 
in the interval < R < L. In addition, for large values 
of p the diffusivity of the particles is almost zero, so we 
neither consider them. Similarly, very small values of Dq 
are disregarded. 

Numerical simulations of eq. show that a statisti- 
cal steady state of the particles distribution is reached. 
We observe this by computing the cluster coefficient of 
the distribution of particles (see below), and observing 
that for long-times it goes to an average constant value, 
which indicates that the spatial distribution of particles 
is statistically stationary (see fig. JIJ). For the statistical 
stationary state we also note that: a) for a wide range of 
values of Do and R the long-time spatial distribution of 
particles is homogeneous if p is in the range < p < I; 
b) for p > 1 the particles distribute non-uniformly for all 



the values of Dq and R considered. Increasing the value 
of p upon the critical value p c — 1 the unhomogeneous 
distribution of particles turns into clusters. However, if 
we keep on increasing p (as already suggested in the for- 
mer paragraph) the effective diffusivity, Di, is almost zero 
and the particles remain almost without moving, so that 
the time scales to observe a nonhomogenous distribution 
become very large. The aggregation proceeds as follows: 
if a number particles are close to each other their ef- 
fective diffusivity diminishes and they do not get much 
appart forming a small cluster. As other particles, which 
are moving randomly in the system, get near this cluster, 
their diffusivity is reduced and they stay in the aggregate. 

No cluster dynamics is observed in the system. The 
groups of particles remain statistically fixed in the space, 
not changing notably their position, since most of the 
particles in any of them have a very small diffusivity. 
Only the particles in the outer parts have a non-negligible 
diffusion coeficient and slightly move, never getting far of 
the cluster and only very rarely leaving the cluster and 
joining to another one. In principle, it could be surpris- 
ing that large values of Dq do not avoid any possibility of 
aggregation in the system. In fact, in the numerical sim- 
ulations one observes that larger values of Dq favour en- 
counters among particles and the clustered steady state 
is reached in a shorter time. Thus, the value of Dq in 
the model only changes the time scale of the system. As 
already mentioned, increasing values of p homogenises 
the spatial distribution of particles, i.e. they are more 
sparsed in the system, and the number of groups in- 
creases. When clusters are formed in the system, their 
typical size strongly depends on the values of the param- 
eters R, p and the initial number of particles. Also, it 
is always observed that the clusters do not periodically 
distribute in space 8] which is what is expected for real- 
istics models of swarming since periodicity is rarely seen 
in these 0. In fig. J5J we show distinct spatial structures 
observed for different values of the parameters. A single 
cluster and several clusters are shown there. 

Clustering is quantified by using an entropy-like mea- 
sure Sm — _ EtiT' n ¥' where M is the num- 
ber of boxes in which we divide the system, and mj 
is the number of particles inside box i. One has that 
< Sm < hi M, so that the minimum value, Sm = 0, is 
obtained when all the particles are in just one of the 
boxes, and the maximum one, InM, is reached when 
nii — N/M for all i, i.e., Sm decreases when there is 
more aggregation. The clustering coefficient is defined as 
Cm = exp(< Sm >t)/M, where < . > t denotes a tempo- 
ral average in steady state conditions, so that when there 
is no clustering Cm ~ 1 and small values of Cm indicate 
a non-uniform distribution of particles. In fig. Q we 
plot Cm versus p for different values of Dq, and always 
a fixed R = 0.05 (similar plots, not shown, are obtained 
for other R's). The transition at p c = 1 is clearly ob- 
served for the different values of Dq. In addition, for a 
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FIG. 1: We see in this figure the approach to the steady 
state. We plot the clustering coefficient, which gives an idea 
of the spatial distribution of particles (see text), vs time for 
different values of Do and the same value of p = 2. From top 
to bottom: D = 10" 5 , D = 1CT 4 , D = 10" 3 . The smaller 
Do the larger the time needed to reach the steady state. 
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FIG. 3: Clustering coefficient versus the parameter p. Each 
plot corresponds to a fixed value of Do- The circles corre- 
spond to Do = 10 -5 , squares to Do = 10 -4 , and diamonds 
to Do = 10~ 3 . The rest of the parameters are R = 0.05 and 
po = 500. The system is divided in a number of M — 144 
boxes and every point in this figure is calculated from a spa- 
tial distribution of particles at time t = 10 7 . 



should be interpreted in the Ito Calculus. This is so 
because it represents the continuous time limit of a dis- 
crete population dynamics model of nonoverlapping gen- 
erations. A clear discussion on this can be consulted in 
the first of the references in |10| . 

Then the mean-field number density equation for 
eq. J5J (particularly appropriate for not too small densi- 
ties) is given by |1TJ| 

^M = Iv 2 (^(x,t)Vx,t)), (3) 

where p(x, t) is the number density of particles in the 
continuum space-time. Taking into account that in the 
continuum limit Nr(i) — -/Vr(x) = J x _ r | <i? ,drp(r, t), the 
final evolution equation for the density in our model is 



FIG. 2: Long-time (10 7 time steps) spatial distribution of 
particles. The parameters are R = 0.05, initial density po = 
500 for all the plots, and a) D = 10" 4 , p = 1.5; b) D 
If) \ p = 2.0; c) Do = 10~ 5 , p = 3.5; d) Do = 10~ 3 , p = 7.5. 
The system size is L — 1. 

fixed Dq one can also see that as p increases (for values 
in the clustering phase, p > p c ) Cm also increases. This 
indicates that the number of clusters in the system gets 
larger as p increases, as can be observed in fig. |2t>) and 
c). 

In order to obtain a continuum description of the model 
to try to understand better its properties, we first take 
the limit At — > in eq. Q and get the Langevin equation 



dx i (i)=/9(x i (t),*)dW i (t), i = l,...,N, (2) 

where dWj(i) is a Wiener proccess and /3(xj(i),t) = 
y/2Do/Nn(i)P. It is very important to note that eq. @ 



dp(x, t) 

at 



A>V 2 



Ix-rKiJ 



drp(r, t) 



(4) 



with initial condition p(x, 0) = N/L 2 . A proper deriva- 
tion of eq. J3J from the interacting particle dynamics, 
eq. gives rise to a multiplicative noise term that 
has been averaged out in eq. (@J. As already mentioned, 
this is a good approximation for not too small densities. 
Moreover, fluctuations in the density equation, which re- 
flect the discrete nature of the particles, seem to have an 
irrelevant role in the pattern formation instability to be 
discussed below (see also 0). 

Particle number is conserved in eq. and also it is 
implicitly assumed that when the density at a point is 
zero within a range R it remains zero (any possible sin- 
gularity is avoided in eq. (J3J). In this density equation it 
is clear that the system is described by a nonlinear diffu- 
sion equation where the effective diffusivity at any point 
is decreased by the total density in a neighbourhood of 
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the point. With respect to other existent nonlinear dif- 
fusion models, e.g. for insect swarming Q, our model 
presents the following crucial differences: the finite range 
of interaction R, and that, as has already been men- 
tioned, in those classical models the motility is generally 
proportional to the local density, at least to the author 
best knowledge. The following non-dimcnsionalization 
s = Dot/R 2 , u = x/R and p = R 2 p transforms eq.(@J 
into 



dp{u,s) 2 
ds 



p(u,s) 



.(/|u-v|<l dv P( V ^) 



(5) 



from which it is clear that D only renormalizes the time- 
scale of the system. 

To study the clustering, let us make a linear sta- 
bility analysis of the stationary homogeneous solution, 
pa = N/L 2 , of eq. @ . We write p(x, t) = p a + etp(x,t) 
where e is a small parameter, and ip(x,t) the space-time 
dependent perturbation, and substitute it in eq. Q). To 
first order in e, ip evolves as [llj : 



9tV(x, t) = 



D 



DoP 



(PottR 2 )p 

ip(r, t) dr 



(6) 



Ix-rKi? 



Then taking an harmonic perturbation ip{x, t) — 
exp(Ai + ik • x) one arrives at the following dispersion 
relation 



X(K) 



D K 2 2pJi(KR) 



( P0 7rR 2 )r 



KR 



I) 



(7) 



being K the modulus of k, and J\ the first-order 
Bessel function. (Non-dimensionalising eq. J7J, or rather 
performing the stability analysis directly to the non- 
dimensional expression eq. (|3J), one has that X(K) — 
(K 2 /{p 7r)P){^j^ - 1), where K = KR and A = 
XR 2 /Do.) The uniform density is unstable, giving rise 
to aggregation, if A is positive. For K = one has that 
A = for all p (and any Do, po and R). This implies that 
the homogeneous density is neutrally stable to uniform 
perturbations which is due to the conservation of the to- 
tal number of particles. As shown in fig. (0} the growth 
rate A (always a real number) is positive for p > 1, while 
the system is stable for p < 1. This is in agreement 
with the numerical results found for the discrete model. 
Moreover, the instability for values of p larger than 1 is 
in a band of wave vectors within the range < K < K + , 
i.e., it is of type II in the classification of 12]. This 
kind of instability appears typically for systems with a 
conservation law, and are characterised by the fact that 
the growth rate, A, vanishes at K=0, and represents the 
onset of aggregation and formation of groups 9]. The 
phenomenology of aggregation already observed in the 




FIG. 4: Dispersion relation, eq. 0, for different values of p. 
In the vertical we plot the (non-dimensional) growth exponent 
A multiplied by (po 7r ) p , and in the horizontal axis K = KR. 
p = 1 is the critical value for the instability. 




FIG. 5: Density distribution at large time. Do = 10 , R = 
0.1, po = 76.23 and p = 2. Lighter colour indicates larger 
density. 



particle model has its perfect counterpart in the disper- 
sion relation eq. (0, being the values of Do, po aud R 
irrelevant for the onset of clustering. Therefore, the den- 
sity equation perfectly explains the aggregation of the 
particles as a deterministic instability of type II. The 
nonlinearities of the model saturate the exponential lin- 
ear growth for p > 1 given by the dispersion relation. It is 
worth mentioning that if an exponentially decreasing de- 
pendence fo the diffusivity would have been considered, 
D{i) oc exp — (pNn(i)), a similar relation dispersion is 
obtained but with a control parameter p dependent on R 
and po, which is also rather realistic. 

We have also performed a numerical simulation of the 
density equation eq. I0J starting from a random initial 
distribution and using a variable time step in order to 
avoid any numerical instability. A long-time non-uniform 
non-periodic pattern is plotted in fig. (0 for D = 10~ 3 , 
R = 0.1, po = 76.23 and p — 2. As expected, but not 
shown, if p < la uniform distribution of the field is 
obtained for these and other parameter values. 

In summary, a new mechanism for clustering of self- 
propelled particles has been presented. It just considers 
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that the motility of any particle decreases with the crowd- 
ing of its surroundings. The model has been studied and 
its number density equation has been obtained, which 
turns out to be a particular form of a nonlinear diffusion 
equation. The relation dispersion indicating the condi- 
tions for clustering were calculated, so that the clustering 
at the particle level can be understood as a deterministic 
instability of the density equation. In addition, the role 
of the parameters of the model becomes clearer at the 
density description level. Regarding a biological applica- 
tion, motile bacteria or plankton particles propel them- 
self, so that they are typically modelled as self-propelled 
particles, and maybe a minimal mechanism like the one 
discussed in this work can explain some of the yet un- 
known causes for their clustering. Ongoin g re search on 
the nature of the transition to clustering P. ll3l| , on other 
possible functional shapes of the motility of the particles, 
and on the influence of a variable interaction range are 
in consideration for the next future. 
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